Dane zawierają informacje o warunkach hodowli i połowch śledzi, zmienną decyzyjną jest długość śledzia. Jej rozkład zbliżony jest do rozkładu normalnego. Na etapie wizualnej eksploracji danych zaobserwowano spadek długości śledzia w czasie. Analiza korelacji nie doprowadziła do stwierdzenia zależności między atrybutami a zmienną decyzyjną, naomiast wykazała wysoką współzależność zagęszczenia Calanus finmarchicus gat. 1 oraz zagęszczenia Calanus helgolandicus gat. 1. W celu predykcji rozmiaru śledzia zastosowano 3 modele uczenia maszynowego: regresję liniową, metodę wektorów nośnych (SVM) oraz eXtreme gradient boosting (XGB). Analiza ważności atrybutów doprowadziła do wniosku, że najbardziej istotną przyczyną spadku długości śledzia była temperatura przy powierzchni wody.
W celu zapewnienia powtarzalności wyników dla uruchamianego wielokrotnie kodu należy ustawić tzw. ziarno. Użyłam ziarna równego 64.
set.seed(64)
Brakujące dane w pliku zostały oznaczone znakiem zapytania. Wczytanie danych z parametrem na pozwala na zastąpienie ich wyrażeniem NA rozpoznawanym przez język R jako brakujące dane.
dataset <- read.csv("sledzie.csv", na=c('?'))
dataset['X'] <- NULL
Kolumna X zawierająca numer próbki nie jest użyteczna w analizie, dlatego została usunięta.
Każdy wiersz danych zawiera jeden atrybut decyzyjny, którym jest długość złowionego sledzia podana w centymetrach length oraz 14 atrybutów predycyjnych, które mozna podzielić na trzy główne grupy:
Wszystkie atrybuty z wyjątkiem miesiąca połowu (xmonth) są atrybutami numerycznymi. Miesiąc połowu jest atrybutem kategorycznym. Dane składają się z 52 582 wierszy.
Poniżej przedstawiono podstawowe statystyki dla każdego z atrybutów:
basic_stats <- simplify2array(summary(dataset))
basic_stats %>%
kable() %>%
kable_styling() %>%
scroll_box(width = "900px")
| length | cfin1 | cfin2 | chel1 | chel2 | lcop1 | lcop2 | fbar | recr | cumf | totaln | sst | sal | xmonth | nao | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Min. :19.0 | Min. : 0.0000 | Min. : 0.0000 | Min. : 0.000 | Min. : 5.238 | Min. : 0.3074 | Min. : 7.849 | Min. :0.0680 | Min. : 140515 | Min. :0.06833 | Min. : 144137 | Min. :12.77 | Min. :35.40 | Min. : 1.000 | Min. :-4.89000 | |
| 1st Qu.:24.0 | 1st Qu.: 0.0000 | 1st Qu.: 0.2778 | 1st Qu.: 2.469 | 1st Qu.:13.427 | 1st Qu.: 2.5479 | 1st Qu.:17.808 | 1st Qu.:0.2270 | 1st Qu.: 360061 | 1st Qu.:0.14809 | 1st Qu.: 306068 | 1st Qu.:13.60 | 1st Qu.:35.51 | 1st Qu.: 5.000 | 1st Qu.:-1.89000 | |
| Median :25.5 | Median : 0.1111 | Median : 0.7012 | Median : 5.750 | Median :21.673 | Median : 7.0000 | Median :24.859 | Median :0.3320 | Median : 421391 | Median :0.23191 | Median : 539558 | Median :13.86 | Median :35.51 | Median : 8.000 | Median : 0.20000 | |
| Mean :25.3 | Mean : 0.4458 | Mean : 2.0248 | Mean :10.006 | Mean :21.221 | Mean : 12.8108 | Mean :28.419 | Mean :0.3304 | Mean : 520367 | Mean :0.22981 | Mean : 514973 | Mean :13.87 | Mean :35.51 | Mean : 7.258 | Mean :-0.09236 | |
| 3rd Qu.:26.5 | 3rd Qu.: 0.3333 | 3rd Qu.: 1.7936 | 3rd Qu.:11.500 | 3rd Qu.:27.193 | 3rd Qu.: 21.2315 | 3rd Qu.:37.232 | 3rd Qu.:0.4560 | 3rd Qu.: 724151 | 3rd Qu.:0.29803 | 3rd Qu.: 730351 | 3rd Qu.:14.16 | 3rd Qu.:35.52 | 3rd Qu.: 9.000 | 3rd Qu.: 1.63000 | |
| Max. :32.5 | Max. :37.6667 | Max. :19.3958 | Max. :75.000 | Max. :57.706 | Max. :115.5833 | Max. :68.736 | Max. :0.8490 | Max. :1565890 | Max. :0.39801 | Max. :1015595 | Max. :14.73 | Max. :35.61 | Max. :12.000 | Max. : 5.08000 | |
| NA | NA’s :1581 | NA’s :1536 | NA’s :1555 | NA’s :1556 | NA’s :1653 | NA’s :1591 | NA | NA | NA | NA | NA’s :1584 | NA | NA | NA |
Pakietem umozliwiającym wygodne przetwarzanie brakujących danych jest imputeTS. Z jego pomocą możemy obejrzeć statyski dotyczące brakujących danych:
stats <- sapply(dataset, statsNA, printOnly = FALSE)
stats %>%
kable() %>%
kable_styling() %>%
scroll_box(width = "900px")
| length | cfin1 | cfin2 | chel1 | chel2 | lcop1 | lcop2 | fbar | recr | cumf | totaln | sst | sal | xmonth | nao | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| lengthTimeSeries | 52582 | 52582 | 52582 | 52582 | 52582 | 52582 | 52582 | 52582 | 52582 | 52582 | 52582 | 52582 | 52582 | 52582 | 52582 |
| numberNAs | 0 | 1581 | 1536 | 1555 | 1556 | 1653 | 1591 | 0 | 0 | 0 | 0 | 1584 | 0 | 0 | 0 |
| percentageNAs | 0% | 3.01% | 2.92% | 2.96% | 2.96% | 3.14% | 3.03% | 0% | 0% | 0% | 0% | 3.01% | 0% | 0% | 0% |
| naGapLongest | NA | 3 | 3 | 3 | 3 | 2 | 3 | NA | NA | NA | NA | 3 | NA | NA | NA |
| naGapMostFrequent | 52582 | 1 | 1 | 1 | 1 | 1 | 1 | 52582 | 52582 | 52582 | 52582 | 1 | 52582 | 52582 | 52582 |
| naGapMostOverallNAs | 52582 | 1 | 1 | 1 | 1 | 1 | 1 | 52582 | 52582 | 52582 | 52582 | 1 | 52582 | 52582 | 52582 |
Procent brakujących danych nie jest wysoki - liczba pustych wartości nie przekracza 3% w każdej z kolumn??? Sprawdzić. Mimo, że nie jest to duży odsetek, brakujące dane mogą negatywnie wpływać na zdolności predykcyjne niektórych modeli uczenia maszynowego.
Istnieje wiele sposobów na uporanie się z brakującymi danymi. Wiersze, które zawierają takie dane można usunąć, wypełnić zawartością bądź zostawić w formie niezmienionej. Najczęstszą strategią, zapobiegającą utracie danych i umożliwiającą przeprowadzenie procesu uczenia maszynowego jest wypełnienie braków danymi syntetycznymi. W tym celu dla kolumn, które zawierają wartości puste wykorzystałam interpolację wartości.
dataset$cfin1 <- na.interpolation(dataset$cfin1)
dataset$cfin2 <- na.interpolation(dataset$cfin2)
dataset$chel1 <- na.interpolation(dataset$chel1)
dataset$chel2 <- na.interpolation(dataset$chel2)
dataset$lcop1 <- na.interpolation(dataset$lcop1)
dataset$lcop2 <- na.interpolation(dataset$lcop2)
dataset$sst <- na.interpolation(dataset$sst)
W tej sekcji przedstawiono rozkład wartości dla poszczególnych atrybutów.
p_length <- ggplot(dataset, aes(length)) +
geom_histogram(binwidth = 0.5) +
scale_x_discrete(name="Długość [cm]", limits= seq(min(dataset$length), max(dataset$length), by=1) ) +
ylab("Liczba obserwacji") +
ggtitle("Długość śledzia") +
theme(plot.title = element_text(hjust = 0.5))
Rozkład wartości długości śledzia jest zbliżony do rozkładu normalnego o średniej ok. 25.5 cm.
p_cfin1 <- ggplot(dataset, aes(cfin1)) + geom_histogram(binwidth = 1.0) + xlab("Zagęszczenie planktonu") +
ylab("Liczba obserwacji") + ggtitle("Calanus finmarchicus gat. 1") + theme(plot.title = element_text(hjust = 0.5))
p_cfin2 <- ggplot(dataset, aes(cfin2)) + geom_histogram(binwidth = 1.0) + xlab("Zagęszczenie planktonu") +
ylab("Liczba obserwacji") + ggtitle("Calanus finmarchicus gat. 2") + theme(plot.title = element_text(hjust = 0.5))
grid.arrange(p_cfin1, p_cfin2, ncol = 2, nrow = 1)
p_chel1 <- ggplot(dataset, aes(chel1)) + geom_histogram(binwidth = 1.0) + xlab("Zagęszczenie planktonu") +
ylab("Liczba obserwacji") + ggtitle("Calanus helgolandicus gat. 1") + theme(plot.title = element_text(hjust = 0.5))
p_chel2 <- ggplot(dataset, aes(chel2)) + geom_histogram(binwidth = 1.0) + xlab("Zagęszczenie planktonu") +
ylab("Liczba obserwacji") + ggtitle("Calanus helgolandicus gat. 2") + theme(plot.title = element_text(hjust = 0.5))
grid.arrange(p_chel1, p_chel2, ncol = 2, nrow = 1)
p_lcop1 <- ggplot(dataset, aes(x=lcop1)) + geom_histogram(binwidth = 1.0) + xlab("Zagęszczenie planktonu") +
ylab("Liczba obserwacji") + ggtitle("Widłonogi gat. 1") + theme(plot.title = element_text(hjust = 0.5))
p_lcop2 <- ggplot(dataset, aes(x=lcop2)) + geom_histogram(binwidth = 1.0) + xlab("Zagęszczenie planktonu") +
ylab("Liczba obserwacji") + ggtitle("Widłonogi gat. 2") + theme(plot.title = element_text(hjust = 0.5))
grid.arrange(p_lcop1, p_lcop2, ncol = 2, nrow = 1)
Z wykresów wynika, że w większości przypadkóW zagęszczenie planktonu ma rozkład prawostronnie skośny. Dla Calanus finmarchicus gat. 1 oraz Widłonogów gat. 1 obserwujemy występowanie wartości odstających tzw. outlierów. Mogą być one skutkiem błędów popełnionych przy zbieraniu danych i negatywnie wpływać na dalszą analizę. Aby nie tracić zupełnie danych z wierszy, dla których niektóre atrybuty są outlierami, zastąpiono je wartościami NA, a następnie przeprowadzono interpolację.
replace_cfin1 <- which(dataset$cfin1 > 20)
replace_lcop1 <- which(dataset$cfin1 > 90)
for (i in replace_cfin1){
dataset[i, "cfin1"] <- NA
}
for (j in replace_lcop1){
dataset[i, "lcop1"] <- NA
}
dataset$cfin1 <- na.interpolation(dataset$cfin1)
dataset$lcop1 <- na.interpolation(dataset$lcop1)
Po przetworzeniu wartości odstających rozkład Calanus finmarchicus gat. 1 i Widłonogów gat. 1 jest bardziej równomierny i prezentuje się następująco:
p_fbar <- ggplot(dataset, aes(fbar)) + geom_histogram(binwidth = 0.05) + xlab("Ułamek pozostawionego narybku")+
ylab("Liczba obserwacji") + ggtitle("Natężenie połowów w regionie") + theme(plot.title = element_text(hjust = 0.5))
p_recr <- ggplot(dataset, aes(recr)) + geom_histogram(binwidth = 50000.0) + xlab("Liczba śledzi")+
ylab("Liczba obserwacji") + ggtitle("Roczny narybek") + theme(plot.title = element_text(hjust = 0.5))
p_cumf <- ggplot(dataset, aes(cumf)) + geom_histogram(binwidth = 0.02) + xlab("Ułamek pozostawionego narybku")+
ylab("Liczba obserwacji") + ggtitle("Łączne roczne natężenie połowów") + theme(plot.title = element_text(hjust = 0.5))
p_totaln <- ggplot(dataset, aes(totaln)) + geom_histogram(binwidth = 50000.0) + xlab("Liczba śledzi")+
ylab("Liczba obserwacji") + ggtitle("Łączna liczba złowionych ryb") + theme(plot.title = element_text(hjust = 0.5))
grid.arrange(p_fbar, p_recr, p_cumf, p_totaln, ncol = 2, nrow = 2)
Dane dotyczące liczby śledzi w łowiskach mają rozkład dużo bardziej równomierny niż dane dotyczące zagęszczenia planktonu. Natężenia połowów mają rozkład wielomodalny.
p_sst <- ggplot(dataset, aes(sst)) + geom_histogram(binwidth = 0.1) + xlab("Temperatura [°C]")+
ylab("Liczba obserwacji")+ ggtitle("Temperatura przy powierzchni wody")+ theme(plot.title = element_text(hjust=0.5))
p_sal <- ggplot(dataset, aes(sal)) + geom_histogram(binwidth = 0.01) + xlab("Poziom zasolenia [Knudsen ppt]")+
ylab("Liczba obserwacji") + ggtitle("Zasolenie wody") + theme(plot.title=element_text(hjust=0.5))
p_xmonth <- ggplot(dataset, aes(xmonth)) + geom_histogram(binwidth = 1.0) + xlab("Numer miesiąca")+
ylab("Liczba obserwacji") + ggtitle("Miesiąc połowu") + theme(plot.title=element_text(hjust=0.5))
p_nao <- ggplot(dataset, aes(nao)) + geom_histogram(binwidth = 0.2) + xlab("Oscylacja [mb]")+
ylab("Liczba obserwacji") + ggtitle("Oscylacja północnoatlantycka") + theme(plot.title=element_text(hjust=0.5))
grid.arrange(p_sst, p_sal, p_xmonth, p_nao, ncol = 2, nrow = 2)
Najczęściej odnotowaną temperaturą wody było 13.6°C. Może być to związane z faktem, że najwięcej danych zostało zebranych w miesiącach letnich - szczególnie w sierpniu. Zasolenie wody koncentruje się głównie w przedziale 35.5 - 35.55 z niewielkimi odstępstwami.
Poniżej przedstawiono moacierz korelacji dla atrybutów.
corrplot(cor(dataset), tl.col="black")
Żaden z atrybutów nie wykazuje wysokiego współczynnika korelacji w stosunku do atrybutu decyzyjnego length. Najwyższą pozytywną korelację można zaobserwować dla par zagęszczenie widłonogów gat. 1 (lcop1) oraz zagęszczenie Calanus helgolandicus gat. 1 (chel1), zagęszczenie widłonogów gat. 2 (lcop2) oraz zagęszczenie Calanus helgolandicus gat. 2 (chel2) a także łączne roczne natężenie połowów w regionie (cumf) oraz natężenie połowów w regionie (fbar). Negatywną korelację wykazuje para łączna liczba ryb złowionych w ramach połowu (totaln) oraz łączne roczne natężenie połowów w regionie (cumf). Przyjrzyjmy się im bliżej.
p_lcop1_chel1 <- ggplot(dataset, aes(x=lcop1, y=chel1)) + geom_point() + geom_smooth(method=lm) + annotate("text", x = 30, y = 87.25,
label = c(paste("wsp. korelacji =", round(cor(dataset$lcop1, dataset$chel1),2))))
p_lcop2_chel2 <- ggplot(dataset, aes(x=lcop2, y=chel2)) + geom_point() + geom_smooth(method=lm) + annotate("text", x = 25, y = 50,
label = c(paste("wsp. korelacji =", round(cor(dataset$lcop2, dataset$chel2),2))))
p_cumf_fbar <- ggplot(dataset, aes(x=cumf, y=fbar)) + geom_point() + geom_smooth(method=lm) + annotate("text", x = 0.15, y = 0.7,
label = c(paste("wsp. korelacji =", round(cor(dataset$cumf, dataset$fbar),2))))
p_totaln_cumf <- ggplot(dataset, aes(x=totaln, y=cumf)) + geom_point() + geom_smooth(method=lm) + annotate("text", x = 800000, y = 0.34,
label = c(paste("wsp. korelacji =", round(cor(dataset$totaln, dataset$cumf),2))))
grid.arrange(p_lcop1_chel1 , p_lcop2_chel2, p_cumf_fbar, p_totaln_cumf, ncol = 2, nrow = 2, top="Wykresy korelacji dla wybranych par zmiennych")
Istotnie, korelacja wskazanych par zmiennych jest wysoka. Ta informacja zostanie wykorzystana na etapie budowy modeli decyzyjnych.
Poniżej predstawiono wykres zmiany długości śledzia w czasie. W celu lepszej wizualizacji zależności, dane zostały przetworzone za pomocą tzw. średniej kroczącej (ang. moving average). Z wykorzystaniem suwaków znajdujących się pod wykresem można zawęzić przedział, z którego były pobrane próbki.
Na wykresie można dostrzec tendencję spadkową - długość śledzi od pewnego momentu zaczęła wyraźnie maleć.
data <- data.frame(dataset)
data$X <- seq.int(nrow(data))
len_smooth = stats::filter(data$length, rep(1/500, 500), side=2)
plot_ly() %>%
add_lines(x=data$X, y=len_smooth, line=list(dash="solid", width = 2.0, color="green")) %>%
layout(
title = "Zmiana długości śledzia w czasie",
xaxis = list(title= "Numer próbki", rangeslider = list()),
yaxis = list(title = "Długość śledzia [cm]"))
Aby skonstruować dobry model uczenia maszynowego, konieczne jest wcześniejsze uważne przygotowanie danych. Wśród atrybutów znajduje się kolumna xmonth oznaczająca numer miesiąca. Mimo tego, że wartości reprezentowane są jako liczby, jest to atrybut kategoryczny. Dlatego jego wartości zostały zamienione za pomocą one-hot encoding - dodano 12 kolumn, z których dla każdego wiersza tylko jedna przyjmuje wartość 1, a pozostałe - 0.
cols <- colnames(dataset)
dataset2 <- cbind(dataset, dummy(dataset$xmonth, sep = "_"))
cols <- c(cols, c("jan", "feb", "mar", "apr", "may", "jun", "jul", "aug", "sep", "oct", "nov", "dec"))
colnames(dataset2) <- cols
dataset2['xmonth'] <- NULL
head(dataset2) %>%
kable() %>%
kable_styling() %>%
scroll_box(width = "900px")
| length | cfin1 | cfin2 | chel1 | chel2 | lcop1 | lcop2 | fbar | recr | cumf | totaln | sst | sal | nao | jan | feb | mar | apr | may | jun | jul | aug | sep | oct | nov | dec |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 23.0 | 0.02778 | 0.27785 | 2.46875 | 21.43548 | 2.54787 | 26.35881 | 0.356 | 482831 | 0.3059879 | 267380.8 | 14.30693 | 35.51234 | 2.8 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| 22.5 | 0.02778 | 0.27785 | 2.46875 | 21.43548 | 2.54787 | 26.35881 | 0.356 | 482831 | 0.3059879 | 267380.8 | 14.30693 | 35.51234 | 2.8 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| 25.0 | 0.02778 | 0.27785 | 2.46875 | 21.43548 | 2.54787 | 26.35881 | 0.356 | 482831 | 0.3059879 | 267380.8 | 14.30693 | 35.51234 | 2.8 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| 25.5 | 0.02778 | 0.27785 | 2.46875 | 21.43548 | 2.54787 | 26.35881 | 0.356 | 482831 | 0.3059879 | 267380.8 | 14.30693 | 35.51234 | 2.8 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| 24.0 | 0.02778 | 0.27785 | 2.46875 | 21.43548 | 2.54787 | 26.35881 | 0.356 | 482831 | 0.3059879 | 267380.8 | 14.30693 | 35.51234 | 2.8 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| 22.0 | 0.02778 | 0.27785 | 2.46875 | 21.43548 | 2.54787 | 26.35881 | 0.356 | 482831 | 0.3059879 | 267380.8 | 14.30693 | 35.51234 | 2.8 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
Kolejnym krokiem jest podział zbioru na zbiór uczący, testowy i walidacyjny. Przyjęto proporcje odpowiednio: 7:2:1.
trainIndex <- createDataPartition(dataset2$length, p=0.7, list=FALSE)
train <- dataset2[trainIndex,]
test <- dataset2[-trainIndex,]
testIndex <- createDataPartition(test$length, p=0.66, list=FALSE)
validate <- test[-testIndex,]
test <- test[testIndex,]
cat("train set: ", 100*nrow(train)/nrow(dataset2), "%", "\n", "test set: ", 100*nrow(test)/nrow(dataset2), "%", "\n", "validation set: ", 100*nrow(validate)/nrow(dataset2), "%")
## train set: 70.00494 %
## test set: 19.79955 %
## validation set: 10.1955 %
Pierwszą wykorzystaną metodą będzie regresja liniowa. Do treningu zastosowano optymalizację według miary RMSE.
fit_lm <- train(length~., data=train, method="lm", metric="RMSE")
Obliczmy błąd popełniany przez model na zbiorze walidacyjnym w ujęciu RMSE oraz R^2.
predictions_lm <- predict(fit_lm, validate)
rmse_lm <- sqrt(mean((validate$length - predictions_lm)^2))
r2_lm <- cor(validate$length, predictions_lm) ^ 2
cat("RMSE: ", rmse_lm, "R2: ", r2_lm)
## RMSE: 1.362388 R2: 0.3292823
Regresja liniowa jest bardzo prostym modelem uczenia maszynowego o ograniczonych możliwościach, jednak błąd popełniony przez tą metodę nie jest bardzo duży.
Kolejnym regresorem wykorzystanym do predykcji długości śledzia jest metoda wektorów nośnych (ang. support vector machines, SVM). Implementacja tej metody znajduje się w pakiecie kernlab.
fit_svm <- ksvm(length~., train)
predictions_svm <- predict(fit_svm, validate)
rmse_svm <- sqrt(mean((validate$length - predictions_svm)^2))
r2_svm <- cor(validate$length, predictions_svm) ^ 2
cat("RMSE: ", rmse_svm, "R2: ", r2_svm)
## RMSE: 1.152424 R2: 0.5200524
Otrzymane wyniki są lepsze niż w przypadku regresji liniowej wedle miary RMSE, natomiast gorsze wedle R^2.W tym przypadku SVM trening trwał stosunkowo długo.
Ostatnią metodą jest eXtreme Gradient Boosting, w skrócie XGB. Wykorzystano implementację dostępną w pakiecie catboost. Metoda ta jest znana ze swojej wysokiej wydajności i efektywności. Wymaga innego przygotowania danych niż poprzednie metody - wymaga osobnego podania predyktorów i zmiennej wyjściowej, a także nie wymaga zakodowania one-hot zmiennych kategorialnych dzięki wbudowanemu mechanizmowi (konnieczność przekazania nazw zmiennych parametrem cat_features). Optymalizacja jest przeprowadzana domyslnie według miary RMSE.
trainIndex_xgb <- createDataPartition(dataset$length, p=0.7, list=FALSE)
train_xgb <- dataset[trainIndex_xgb,]
test_xgb <- dataset[-trainIndex_xgb,]
testIndex_xgb <- createDataPartition(test_xgb$length, p=0.66, list=FALSE)
validate_xgb <- test_xgb[-testIndex_xgb,]
test_xgb <- test_xgb[testIndex_xgb,]
train_Y <- train_xgb$length
test_Y <- test_xgb$length
validate_Y <- validate_xgb$length
train_xgb$length <- NULL
test_xgb$length <- NULL
validate_xgb$length <- NULL
train_pool <- catboost.load_pool(train_xgb, label=train_Y, cat_features = c("xmonth"))
test_pool <- catboost.load_pool(test_xgb, label=test_Y, cat_features = c("xmonth"))
val_pool <- catboost.load_pool(validate_xgb, label=validate_Y, cat_features = c("xmonth"))
fit_params <- list(iterations = 10000)
fit_catboost <- catboost.train(train_pool, test_pool = val_pool, params=fit_params)
predictions_xgb <- catboost.predict(fit_catboost, val_pool)
rmse_xgb <- sqrt(mean((validate$length - predictions_xgb)^2))
r2_xgb <- cor(validate$length, predictions_xgb) ^ 2
cat("RMSE: ", rmse_xgb, "R2: ", r2_xgb)
## RMSE: 1.591494 R2: 0.1752335
Uzyskane wyniki są najlepsze z dotychczas wykorzystanych metod uczenia maszynowego pod kątem miary R^2, jednak gorsze pod względem miary RMSE. Jednakże błąd R^2 jest prawie 2 razy mniejszy niż dla najlepszego dotychczasowego wyniku wedle tej miary, a RMSE tylko nieznacznie gorsze, dlatego dalszej analizie poddano tę metodę.
Na etapie analizy korelacji wyłoniono parę lcop1 i chel1 jako zmienne charakteryzujące sie najwyższym współczynnikiem korelacj równym 0.96. Spróbujmy usunąć zmienną chel1 i zaobserwować wpływ tego zabiegu na trafność predykcji modelu.
fit_params <- list(iterations = 10000, ignored_features=c(2))
fit_catboost_ch1 <- catboost.train(train_pool, test_pool = val_pool, params=fit_params)
predictions_xgb_ch1 <- catboost.predict(fit_catboost_ch1, val_pool)
rmse_xgb_ch1 <- sqrt(mean((validate$length - predictions_xgb_ch1)^2))
r2_xgb_ch1 <- cor(validate$length, predictions_xgb_ch1) ^ 2
cat("RMSE: ", rmse_xgb_ch1, "R2: ", r2_xgb_ch1)
## RMSE: 1.588176 R2: 0.1766239
Usunięcie chel2 ze zbioru atrybutów nie przyniosło znaczącej różnicy. Oznacza to jednak, że dla tej metody można wykorzystać mniej zmiennych bez utraty jakości predykcji, co jest dobrą informacją.
Na koniec przetestowano zbudowane modele na zbiorze testowym.
test_lm <- predict(fit_lm, test)
test_svm <- predict(fit_svm, test)
test_xgb <- catboost.predict(fit_catboost, test_pool)
test_xgb_ch1 <- catboost.predict(fit_catboost_ch1, test_pool)
rmse_test_lm <- sqrt(mean((test$length - test_lm)^2))
r2_test_lm <- cor(test$length, test_lm) ^ 2
rmse_test_svm <- sqrt(mean((test$length - test_svm)^2))
r2_test_svm <- cor(test$length, test_svm) ^ 2
rmse_test_xgb <- sqrt(mean((test$length - test_xgb)^2))
r2_test_xgb <- cor(test$length, test_xgb) ^ 2
rmse_test_xgb_ch1 <- sqrt(mean((test$length - test_xgb_ch1)^2))
r2_test_xgb_ch1 <- cor(test$length, test_xgb_ch1) ^ 2
test_res <- data.frame("metoda" = c("regresja liniowa","SVM","XGB", "XGB bez chel1"),
"RMSE" = c(rmse_test_lm, rmse_test_svm, rmse_test_xgb, rmse_test_xgb_ch1),
"R2" = c(r2_test_lm, r2_test_svm, r2_test_xgb, r2_test_xgb_ch1))
test_res %>%
kable() %>%
kable_styling()
| metoda | RMSE | R2 |
|---|---|---|
| regresja liniowa | 1.352261 | 0.3315246 |
| SVM | 1.153885 | 0.5133021 |
| XGB | 1.513417 | 0.2246637 |
| XGB bez chel1 | 1.510245 | 0.2259220 |
Wyniki na zbiorze testowym wykazują podobne względne różnice w jakości predykcji jak dla zbioru walidacyjnego.
Implementacja catboost pozwala na analizę ważności cech atrybutów. Przedstawia listę atrybutów wraz z ich procentową kontrybucją do ostatecznego wyniku predykcji.
catboost.get_feature_importance(fit_catboost)
## cfin1 cfin2 chel1 chel2 lcop1 lcop2 fbar
## 7.711926 7.286647 4.745681 5.103176 4.424734 5.706697 8.235664
## recr cumf totaln sst sal xmonth nao
## 6.253652 5.489162 6.940089 17.395016 5.807484 6.963496 7.936575
Można zaobserwować, że najbardziej istotnym atrybutem okazała się temperatura przy powierzchni wody, która stanowi około 15% wkładu do predykcji.